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ABSTRACT 

We have performed high-resolution 3D simulations of relativistic jets with 
beam flow Lorentz factors up to 7, a spatial resolution of 8 cells per beam radius, 
and for up to 75 normalized time units to study the morphology and dynamics 
of 3D relativistic jets. Our simulations show that the coherent fast backfiows 
found in axisymmetric models are not present in 3D models. We further find 
that when the jet is exposed to non-axisymmetric perturbations, (i) it does not 
display the strong perturbations found for 3D classical hydrodynamic and MHD 
jets (at least during the period of time covered by our simulations), and (ii) it 
does propagate according to the ID estimate. Small 3D effects in the relativistic 
beam give rise to a lumpy distribution of apparent speeds like that observed 
in M87. The beam is surrounded by a boundary layer of high specific internal 
energy. The properties of this layer are briefly discussed. 

Subject headings: galaxies: jets — hydrodynamics — methods: numerical — 
relativity 
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1. Introduction 

Since several years tlie dynamical and morphological properties of axisymmetric 
relativistic jets are investigated by means of relativistic hydrodynamic simulations (van 
Putten 1993; Duncan & Hughes 1994; Marti et al. 1994, 1995, 1997; Komissarov & Falle 
1998; Rosen et al. 1999). In addition, relativistic MHD simulations have been performed in 
2D (Koide, Nishikawa & Muttel 1996; Koide 1997) and 3D (Nishikawa et al. 1997, 1998). In 
their 3D simulations Nishikawa et ai. have studied mildly relativistic jets (Lorentz factor 
4.56) propagating both along and obliquely to an ambient magnetic field. In this Letter 
we report on high-resolution 3D simulations of relativistic jets with the largest beam flow 
Lorentz factor performed up to now (7.09), the largest resolution (8 cells per beam radius), 
and covering the longest time evolution (75 normalized time units; a normalized time unit is 
defined as the time needed for the jet to cross a unit length; see Massaglia, Bodo & Ferrari 
1996). 

The calculations have been performed with the high-resolution 3D relativistic 
hydrodynamics code GENESIS (Aloy et ai. 1999), which is an upgraded version of the code 
developed by Marti, Miiller & Ibanez (1994) and Marti et al. (1995). GENESIS integrates 
the 3D relativistic hydrodynamic equations in conservation form in Cartesian coordinates 
including an additional conservation equation for the density of beam material. The 
computations were performed on a Cartesian domain (X,Y,Z) of size 15Ri, x 15-Rf, x 75i?6 
(120 X 120 X 600 computational cells), where Rf, is the beam radius. The jet is injected at 
z = along the positive z-axis through a circular nozzle defined by x"^ + y"^ < Rf. Beam 
material is injected with a beam mass fraction / = 1, and the computational domain is 
initially filled with an external medium (/ = 0). 

We have considered a 3D model corresponding to model C2 of Marti et al. (1997), 
which is characterized by a beam-to-external proper rest-mass density ratio rj = 0.01, a 
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beam Mach number M;, = 6.0, and a beam flow speed v^ = 0.99c (c is the speed of liglit) 
or a beam Lorentz factor Wf, ~ 7.09. An ideal gas equation of state with an adiabatic 
exponent 7 = 5/3 is assumed to describe both the jet matter and the ambient gas. The 
beam is assumed to be in pressure equilibrium with the ambient medium. 

The evolution of the jet was simulated up to T si ISO-R^/c, when the head of the jet is 
about to leave the grid. The mean jet propagation speed Vh ~ 0.5c, while the ID estimate 
of the jet propagation speed (see, e.g., Marti et al. 1997) gives 0.42c, i.e., our simulations are 
still within the ID phase (see Marti, Miiller & Ibanez 1998). 

Non-axisymmetry was imposed by means of a helical velocity perturbation at the 
nozzle given by 

< = (Vbcos ( — ] , vf = (vbsm ( — ] , vl = i;feyl^-C^, (1) 

where ( is the ratio of the toroidal to total velocity and r the perturbation period 



i.e.,T = T/n, n being the number of cycles completed during the whole simulation). The 
i^avelength of the perturbation. A, is 
L is the axial dimension of the grid. 



wavelength of the perturbation. A, is obtained from the expression A = f^r ^ , where 



VhU 



2. Morphology and dynamics of 3D relativistic jets 

We have considered a model with a 1% perturbation in helical velocity {C, = 0.01) and 
n = 50. Figure|l| shows various quantities of the jet in the plane y = at the end of the 
simulation. Two values of the beam mass fraction are marked by white contour levels. 
The beam structure is dominated by the imposed helical pattern with a characteristic 
wavelength of ~ S.ORf, (to be compared with the value A = 3.5-Rf, expected from the 
estimate of A in the previous paragraph) and an amphtude of ~ 0.2/?^. 
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2.1. Cocoon 

The overall jet's morphology is characterizad by the presence of a highly turbulent, 
subsonic, asymmetric cocoon. The pressure distribution outside the beam is nearly 
homogeneous giving rise to a symmetric bow shock (Fig.|l|b). As in the classical case 
(Norman 1996), our relativistic 3D simulation shows less ordered structures in the cocoon. 
The cocoon remains quite thin (~ 2i?b) as long as the jet propagates efficiently. 

The flow field outside the beam shows that the high velocity backfiow is restricted 
to a small region in the vicinity of the hot spot (Fig.|l|e), the largest backfiow velocities 
(~ 0.5c) being significantly smaller than in 2D models. The flow with high Lorentz factor 
found in axisymmetric simulations (see flow patterns in Marti et al. 1996) appears here 
restricted to a thin layer around the beam and possesses sub-relativistic speeds (~ 0.25c). 
The magnitude of the backfiow velocities in the cocoon do not support relativistic beaming. 

2.2. Beam and hot spot 

Within the beam the perturbation pattern is superimposed to the conical shocks at 
about 26 and 50 Rb- The beam does not exhibit the strong perturbations (deflection, 
twisting, flattening or even filamentation) found by other authors (Norman (1996) for 3D 
classical hydrodynamic jets; Hardee (1996) for 3D classical MHD jets). This can be taken as 
a sign of stability, although it can be argued that our simulation is not evolved far enough. 
Obviously, the beam cross section and the internal conical shock structure are correlated 
(bottom panel in Figure 0). Before the ffist recoUimation shock the beam cross section 
shrinks to an effective radius of 0.7/?;,. After this shock and in the rarefaction the beam 
reexpands and stretches due to an elliptical surface mode (e.g., Hardee 1996). Between 
37Rh^z'^50Rb the beam flow is influenced by the second recoUimation shock, which causes 
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a compression of the beam. A triangular mode seems to grow in this region. 

The hehcal pattern propagates along the jet at nearly the beam speed (see animation 
at [http:/ /scry. uv.es/aloy.html/ J ETS/videos/n50p01|) which could yield to superluminal 
components when viewed at appropriate angles. Besides this superluminal pattern, the 
presence of emitting fluid elements moving at different velocities and orientations could lead 
to local variations of the apparent superluminal motion within the jet. This is shown in 
Fig.^, where we have computed the mean (along each line of sight, and for a viewing angle 
of 40 degrees) local apparent speed. The distribution of apparent motions is inhomogeneous 
and resembles that of the observed individual features within knots in M87 (Biretta, Zhou, 
& Owen 1995). 

The jet can be traced continuously up to the hot spot which propagates as a strong 
shock through the ambient medium. Beam material impinges on the hot spot at high 
Lorentz factors. We could not identify a terminal Mach disk in the flow. We find flow 
speeds near (and in) the hot spot much larger than those inferred from the one dimensional 
estimate. This fact was already noticed for 2D models by Komissarov & Falle (1996) and 
suggested by them as a plausible explanation for an excess in hot spot beaming. 

2.3. Beam/cocoon shear layer 

We find a layer of high specific internal energy (Fig. |l]d) surrounding the beam like in 
previous axisymmetric models (Aloy et al. 1999). A comparison with the backfiow velocities 
(Fig. ||e) shows that it is mainly composed of forward moving beam material at a speed 
smaller than the beam speed. The intermediate speed of the layer material is due to shear 
in the beam/cocoon interface, which is also responsible for its high specific internal energy. 
The existence of such a boundary layer has been invoked by several authors (Komissarov 
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1990, Laing 1996) to interpret a number of observational trends in FRI radio sources. 
Swain, Bridle & Baum (1998) have found evidence for these boundary layers in FRIIs 
(3C353). 

The diffusion of vorticity caused by numerical viscosity is responsible for the formation 
of the boundary layer. Although being caused by numerical effects (and not by the physical 
mechanism of turbulent shear) the properties of PPM-based difference schemes are such 
that they can mimic turbulent flow to a certain degree (Porter & Woodward 1994). Hence, 
our calculations represent a first approach to study the development of shear layers in 
relativistic jets and their observable consequences. The structure of both the shear layer 
and the beam core are sketched in Fig|^. The specific internal energy of the gas in the shear 
layer (region with 0.2 < / < 0.8) is typically more than one order of magnitude larger than 
that of the gas in the beam core. The shear layer broadens with distance from 0.2Rf, near 
the nozzle to l.lRh near the head of the jet (Fig.^. 

2.4. Jet propagation efficiency and disruption 

From the head's position at the end of the simulation (T = 140.8) a mean jet advance 
speed of 0.47c is obtained, but the jet's propagation proceeds in two distinct phases: (i) for 
t^lOO the jet propagates roughly at the estimated ID speed (0.42c); (ii) for t^lOO the jet 
accelerates and propagates at a considerably larger speed (0.55c). Comparing with the 3D 
simulation of Norman (1996) we find a similar behaviour: after a short ID phase and before 
the deceleration, the jet transiently accelerates to a propagation speed which is ~ 20% larger 
than the corresponding ID estimate. This result contradicts the one obtained by Nishikawa 
et al. (1997, 1998), who found a propagation speed of only 70% of the corresponding ID 
estimate in a shorter {^ 20 normalized time units) simulation of a denser jet. Although the 
estimate do not take into account the effect of an extra magnetic pressure in the external 



medium opposing to the jet propagation, in the case of Nishikawa et al. simulations this 
magnetic pressure is neghgible in comparison with the beam momentum density. 

Figure^ shows the axial component of the momentum of the beam particles (integrated 
across the beam) along the axis, which decreases by 30% within the first 60 Rfy. Neglecting 
pressure and viscous effects, and assuming stationarity the axial momentum should be 
conserved, and hence the beam flow is decelerating. The momentum loss goes along with 
the growth of the boundary layer whose material is accelerated and heated by viscous 
stresses. Biconical shocks in the beam are responsible for the break in the axial momentum 
profiles at z = 26Rb and z = 50Rb, because when the beam material passes a conical shock 
and enters into the adjacent rarefaction fan, it is accelerated by local pressure gradients. 

How can the jet accelerate while the beam material is decelerating? Although the beam 
material decelerates, its terminal Lorentz factor is still large enough to produce a fast jet 
propagation. On the other hand, in 3D, the beam is prone to strong perturbations which 
can affect the jet's head structure. In particular, a simple structure like a terminal Mach 
shock will probably not survive when significant 3D effects develop. It will be substituted 
by more complex structures in that case, e.g., by a Mach shock which is no longer normal 
to the beam flow and which wobbles around the instantaneous flow direction. Another 
possibility is the generation of oblique shocks near the jet head due to off-axis oscillations of 
the beam. Although difficult to check quantitatively (due to both the lack of an operative 
definition for Mach disk identification and the present resolution of our simulations) both 
possibilities will cause a less efficient deceleration of the beam fiow at least during some 
epochs. At longer time scales the growth of 3D perturbations will cause the beam to spread 
its momentum over a much larger area than that it had initially, which will efficiently 
reduce the jet advance speed. 
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3. Conclusions 

We have presented a first attempt to analyze the morpho-dynamical properties of 3D 
relativistic jets. From our simulations, we can conclude that the coherent fast backfiows 
found in axisymmetric models are not present in 3D models. We have investigated the 
beam's response to non-axisymmetric perturbations to check its stability. During the period 
of time studied by us (t^l50i?b/c), the beam does not display the strong perturbations 
found by other authors in classical jets (Norman 1996, Hardee 1996) and propagates 
according to the ID estimate. Small 3D effects in the relativistic beam give rise to a lumpy 
distribution of apparent speeds like that observed in M87 (Biretta, Zhou & Owen 1995). 
We have also analyzed the properties of the boundary layer present in our model. 

Obviously, our study must be extended to a wider range of models and perturbations. 
In particular, stronger perturbations should be considered to reach the nonlinear regime 
and to identify the acoustic and mixing phases (Bodo 1998) leading to the jet disruption. 
Further investigation also requires the dependence of the shear layer properties on the 
perturbation parameters. Finally, appropriate perturbations can be studied that mimic the 
wiggles observed in specific sources both at pc (0836+710, Lobanov et al. 1998; 0735+178, 
Gomez et ai. 1999) and kpc scales (M87; Biretta, Zhou & Owen 1995). 

ACKNOWLEDGEMENTS 

This work has been supported in part by the Spanish DGES (grants PB97-1164 
and PB97-1432) and the CSIC-Max Planck Gesellschaft agreement. MAA expresses 
his gratitude to the Conselleria d'Educacio i Ciencia de la Generalitat Valenciana for 
a fellowship. The calculations were carried out on two SGI Origin 2000 at the Centre 
Europeu de Paral.lelisme de Barcelona and at the Centre de Informatica de la Universitat 
de Valencia. 



-10- 

REFERENCES 

Aloy, M.A., Ibafiez, J.M^., Marti, J.M^. & Miiller, E. 1999, ApJS, in press 

Biretta, J.A., Zhou, F., & Owen, F.N. 1995, ApJ, 447, 582 

Bodo, G., 1998, in Astrophysical jets: Open problems, eds. S. Massaglia & G. Bodo, Gordon 
and Breach, p. 161 

Duncan, G.C., & Hughes, P.A. 1994, ApJ, 436, L119 

Gomez, J.L., Marscher, A. P., Alberdi, A., & Gabuzda, D.C. 1999, ApJ, in press 

Hardee, P.E. 1996, in Energy Transport in Radio Galaxies and Quasars, eds. P.E. Hardee, 
A.H. Bridle, J. A. Zensus. ASP Conference Series, Vol. 100, p. 273 

Koide, S. 1997, ApJ, 487, 66 

Koide, S., Nishikawa, K.-L, & Muttel, R.L. 1996, ApJ, 463, L71 

Komissarov, S.S. 1990, Sov. Astron. Lett. 16(4), 284 

Komissarov, S.S., & Falle, S.A.E.G. 1996, in Energy Transport in Radio Galaxies and 
Quasars, eds. P.E. Hardee, A.H. Bridle, J. A. Zensus. ASP Conference Series, Vol. 
100, p. 327 

Komissarov, S.S., & Falle, S.A.E.G. 1998, MNRAS, 297, 1087 

Laing, R.A. 1996, in Energy Transport in Radio Galaxies and Quasars, eds. P.E. Hardee, 
A.H. Bridle, J. A. Zensus. ASP Conference Series, Vol. 100, p. 241 

Lobanov, A. P., Krichbaum, T.P., Witzel, A., Kraus, A., Zensus, J. A., Britzen, S., Otterbein, 
K., Hummel, C.A., & Johnston, K. 1998, A&A, 340, L60 



- 11- 

Marti, J.M^., Font, J.A., Ibafiez, J.M^.& Miiller, E. 1996, in Energy Transport in Radio 
Galaxies and Quasars, eds. P.E. Hardee, A.H. Bridle, J. A. Zensus. ASP Conference 
Series, Vol. 100, p. 149 

Marti, J.M^, Miiller, E., Font, J. A., & Ibaiiez, J.M^ 1995, ApJ, 448, L105 

Marti, J.M^., Miiller, E., Font, J.A., Ibaiiez, J.M^. & Marquina, A. 1997, ApJ, 479, 151 

Marti, J. M^, Miiller, E., & Ibaiiez, J. M^ 1994, A&A, 281, L9 

Marti, J.M— ., Miiller, E. & Ibaiiez, J.M— . 1998, in Astrophysical jets: Open problems, eds. 
S. Massaglia & G. Bodo, Gordon and Breach, p. 149 

Massaglia, S., Bodo, G., & Ferrari, A. 1996, A&A, 307, 997 

Nishikawa, K.-L, Koide, S., Sakai, J., Christodoulou, D.M., Sol, H., & Mutel, R.L. 1997, 
ApJ, 483, L45 

Nishikawa, K.-L, Koide, S., Sakai, J., Christodoulou, D.M., Sol, H., & Mutel, R.L. 1998, 
ApJ, 498, 166 

Norman, M.A. 1996, in Energy Transport in Radio Galaxies and Quasars, eds. P.E. Hardee, 
A.H. Bridle, J. A. Zensus. ASP Conference Series, Vol. 100, p. 319 

Porter, D.H. & Woodward, P.R. 1994, ApJS, 93, 309 

Rosen, A., Hughes, P. A., Duncan, G. C, & Hardee, P. E. 1999, AJ, in press 

Swain, M.R., Bridle, A.H., & Baum, S.A. 1998, ApJ, 507, L29 

van Putten, M.H.P.M. 1993, ApJ, 408, L21 



This manuscript was prepared with the AAS L^T[t;X macros v4.0. 



-12- 





Fig. 1. — Rest-mass density, pressure, flow Lorentz factor, specific internal energy and 
backflow velocity distributions (from top to bottom) of the model discussed in the text in 
the plane y = at the end of the simulation. White contour levels appearing in each frame 
correspond to values of / equal to 0.95 (inner contour; representative of the beam) and 
0.05 (representative of the cocoon/shocked external medium interface). The bottom panel 
displays the isosurface of / = 0.95. 



-13- 




0.9 



Fig. 2. — Mean local apparent speed for the jet of Fig. 1 observed at an angle of 40 degrees. 
Arrows show the projected direction and magnitude of the apparent motion the contours 
corresponding to values of 1.0 c, 1.6 c, and 2.2 c, respectively. Averages have been computed 
along the line of sight for each pixel in the image (computational cell) using the emission 
coefficient as a weight. 
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Structure of the beam/cocoon shear layer 
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Fig. 3. — Beam mass fraction (dotted line), flow Lorentz factor (filled line) and specific 
internal energy (dashed line), in arbitrary units, accross the beam [z = 11. 7 Rb). 



-15- 



Beam mean radius and axial momentum 
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Fig. 4. — Dashed line: Axial component of the momentum of the beam particles (integrated 
accross the beam) along the jet axis at the end of the simulation. Solid lines: Mean beam 
radius along the axis for / > 0.2 (top line) and / > 0.8 (bottom line), respectively. Quantities 
are in code units. 



